Working with Spatial Data in terra and sf

sf is a new package and appears to be the leading contender to replace rgdal and sp. sf reads and writes shapefiles and also has a large number of vector transforms available (see resources at the bottom of the page).

In the code below, you'll notice that R orders geographic coordinate values as "longitude, latitude" rather than lat,lon. This is common in spatial programming as we order everything as x,y and typically we treat longitude as an x and latitude as a y when working with geographic data. This also alighns with Easting, Northing.

The rspatial.org website is the best resource I have found for this material but it does not have the best tutorials. There are a lot of other websites out there but be careful as many of them are out of date.

Working with Vector Data

Loading a Shapefile

The code below will read a shapefile into a dataframe.

# Load libraries
library(terra)
library(sf)

# set the working directory so we don't have to specify a full path each 
# time we load/save a file
setwd("C:/Users/Jim/Desktop/GSP 570/Lab 4 GLM")

# load the shapefile into a SpaVector object
TerraVectorObject=vect("test.shp")

# use the standard plot function to show the spatial data
plot(TerraVectorObject)

# get the number of features in the shapefile
NumFeatures=nrow(TerraVectorObject)

# get one column of attribute values 
ColumnOfAttributeValues=TerraVectorObject$measyear

Because a SpatVector object is similar to a dataframe, the dataframe functions such as nrow() and ncol() work with the object.

Saving a Shapefile

writeVector(TheSpatVectorObject, "TestWriteTerra.shp")

Loading a CSV with Points

TheSpatVectorObject=vect(DataTableOfPoints,geom=c("lon", "lat"),crs=CRSDef)
plot(TheSpatVectorObject)

Saving a SpatVector object to a CSV

If you have a SpatVector object, you'll need to convert it to a data frame and then write it to a CSV file:

TheDataFrame=as.data.frame(TerraVectorObject, geom="XY") 
write.csv(TheDataFrame,"C:/Temp/testy.csv")

Accessing the Spatial Data

The SpatVector object handles points, polylines, and polygons. We'll mostly be using points. We can access the spatial data using the 'geom()' function. This will return a dataframe that contains the x (lon or easting) and y (lat or northing) values in columns 3 and 4. Then, we can access each feature's coordinate by indexing into the dataframe as below.

TheGeom=geom(TerraVectorObject) # provides a matrix of values

# get the coorindate for the second feature in the file
Lon=TheGeom[2,3] # column 3 are the lons
Lat=TheGeom[2,4] # column 4 contains the lats

ProjectingVector Objects

Projecting SpatVector and SpatRaster objects is easy once you have the proj string to define the projection.

ProjectedRaster=project(CroppedRaster, "+proj=utm +zone=10")
plot(ProjectedRaster)

Coordinate Reference Systems

You can use Proj strings to create Coordinate Reference Systems (crs) (a.k.a. Spatial Reference Systems). The example below creates an unproject (geographic) crs and then create a vector dataset that uses that crs.

CRSDef <- "+proj=longlat +datum=WGS84" # proj string

pts <- vect(lonlat, crs=crdref)

Other proj strings include:

+proj=utm +zone=10 # UTM Zone 10 North (North is the default)
+proj=utm +zone=59 +south # UTM Zone 59 South

# Robinson projection method. Try changing the longitude of origin
+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m # robinson in meters

The state plane CRSs have different projection methods and parameters. This information is available online at a number of websites (see below). As an example, State Plane California uses the Lambert Conic Conformal project method. Which has the following parameters (maptools.org).

+proj=lcc   
+lat_1=Latitude of first standard parallel             
+lat_2=Latitude of second standard parallel               
+lat_0=Latitude of false origin (i.e. latitude of origin)                 
+lon_0=Longitude of false origin (i.e. longitude of origin)
+x_0=False Origin Easting                
+y_0=False Origin Northing

The proj string for NAD 83 California State Plane Zone I meters would then be (EPSG.io):

+proj=lcc +lat_1=41.6666666666667 +lat_2=40 +lat_0=39.3333333333333 
+lon_0=-122 +x_0=6561666.667 +y_0=1640416.667

You may see the "+no_defs" parameter in some examples. This is an old parameter that prevented proj from reading a defualts file. This feature no longer exists in proj so "no_defs" it not needed.

You also may see "+towgs84=0,0,0" or something similar. This defines a shift between datums but if the values are 0 then it is not needed.

Additional Resources

Simple Features for R - website for SF with great quick reference.

Spatial Data Science - Vector Data - great detailed explaination of tarra.